Thierry Phénix & Ladislas Nalborczyk
17 mars 2018
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
Interpréter ce genre de données
Construire ce type de modèle
\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i}&= \alpha + \beta x_{i} \\ \alpha &\sim \mathrm{Normal}(60,10) \\ \beta &\sim \mathrm{Normal}(0,10) \\ \sigma &\sim \mathrm{HalfCauchy}(0,1) \end{align} \]
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton
d <- read.csv("Data/Howell1.csv")
str(d)
'data.frame': 544 obs. of 5 variables:
$ X : int 1 2 3 4 5 6 7 8 9 10 ...
$ height: num 152 140 137 157 145 ...
$ weight: num 47.8 36.5 31.9 53 41.3 ...
$ age : num 63 63 65 41 51 35 32 27 19 54 ...
$ male : int 1 0 0 1 0 1 0 1 0 1 ...
d2 <- d %>% filter(age >= 18)
head(d2)
X height weight age male
1 1 151.765 47.82561 63 1
2 2 139.700 36.48581 63 0
3 3 136.525 31.86484 65 0
4 4 156.845 53.04191 41 1
5 5 145.415 41.27687 51 0
6 6 163.830 62.99259 35 1
d2 <- d %>% filter(age >= 18)
head(d2)
X height weight age male
1 1 151.765 47.82561 63 1
2 2 139.700 36.48581 63 0
3 3 136.525 31.86484 65 0
4 4 156.845 53.04191 41 1
5 5 145.415 41.27687 51 0
6 6 163.830 62.99259 35 1
set.seed(1234)
tibble(value = rnorm(1000, 155, 10)) %>%
ggplot(aes(x = value) ) + xlim(130, 180) +
geom_histogram(bins = 20, col = "white") +
theme_bw(base_size = 30)
La likelihhod caractérisant les hauteurs \( ~h_{i}~ \) peut-être représentée par le modèle suivant:
\[ h_{i} \sim \mathrm{Normal}(\mu,\sigma) \]
Elle admet deux paramètres \( ~\mu~ \) et \( ~\sigma \).
On cherche à calculer LA distribution postérieure qui décrit le mieux la répartition des tailles \( ~\{h_i\} \).
\[ \color{purple}{p(\mu,\sigma~|~\{h_i\})} \]
\( \longrightarrow \) On va donc explorer toutes les combinaisons possibles de \( \mu \) et \( \sigma \) et les classer par leurs probabilités respectives.
Notre but est de décrire la distribution a posteriori, qui sera donc une distribution de distributions.
Le calcul de la distribution postérieure des hauteurs est le produit de la fonction de vraisemblance et de la distribution a priori:
\[ \color{steelblue}{p(\mu, \sigma) = p(\mu)~p(\sigma)} \]
\[ \color{orangered}{f(h_i~|~\mu, \sigma)=\frac{1}{\sigma\sqrt{2\pi}} \exp\bigg(-\frac{1}{2}\big(\frac{h_i-\mu~}{\sigma}\big)^{2}\bigg)} \]
\( \color{steelblue}{\mu \sim \mathrm{Normal}(178,20)} \)
tibble(
x = seq(from = 100, to = 250, by = .1),
y = dnorm(x, mean = 178, sd = 20) ) %>%
ggplot(aes(x, y)) +
geom_line( col = "steelblue", lwd = 2) +
xlab(expression(mu)) + ylab("density") +
theme_bw(base_size = 30)
\( \color{steelblue}{\sigma \sim \mathrm{Uniform}(0,50)} \)
tibble(
x = seq(from = -10, to = 60, by = .1),
y = dunif(x, 0, 50)) %>%
ggplot(aes(x, y)) +
geom_line( col = "steelblue", lwd = 2) +
xlab(expression(sigma)) + ylab("density") +
theme_bw(base_size = 30)
Calcul 'force brute' méthode grille
prior_p <- function(x, y) {
dnorm(x, mean = 178, sd = 20)*dunif(y, 0, 50)
}
n <- 200
mu <- seq(from = 100, to = 250, length.out = n)
sigma <- seq(from = -10, to = 60, length.out = n)
prior <- outer(mu, sigma, prior_p)
prior <- prior / sum(prior)
persp3Drgl(mu, sigma, prior, smooth = FALSE, lighting = TRUE, new = TRUE)
You must enable Javascript to view this page properly.
set.seed(1234)
n <- 1e4
tibble(
sample_mu = rnorm(n, mean = 178, sd = 20),
sample_sigma = runif(n, min = 0, max = 50)) %>%
mutate(x = rnorm(n, mean = sample_mu, sd = sample_sigma)) %>%
ggplot(aes(x = x)) +
geom_density(colour = "steelblue", size = 0) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = expression(paste("Prior predictive distribution for ", italic(h[i]))),
x= NULL) +
theme_bw(base_size = 30)
\[ p(h_i = 162.8648~|~\mu, \sigma) \] \[ = 0~~~ \forall \mu, \sigma \]
\[ f(h_i = 162.8648]~|~\mu = 151.23, \sigma = 23.42) \] \[ = 0.01505675 \]
où la \( 34^{ème} \) mesure de hauteur vaut \( 162.8648 \)
\[ \begin{align*} \color{purple}{p(\mu, \sigma~|~\{h_i\})} & = \frac{\color{steelblue}{\mathrm{Normal}(\mu|178,20)\mathrm{Uniform}(\sigma|0,50)}~\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i}|\mu, \sigma)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i}|\mu,\sigma)\mathrm{Normal}(\mu|178,20)\mathrm{Uniform}(\sigma|0,50)d\mu d\sigma}} \\ \\ & \propto \color{steelblue}{\mathrm{Normal}(\mu|178,20)\mathrm{Uniform}(\sigma|0,50)}~\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i}|\mu, \sigma)} \end{align*} \]
Il s'agit de la même formule vue lors des cours 1 et 2, mais cette fois en considérant qu'il existe plusieurs observations de taille (\( h_{i} \)), et deux paramètres à estimer \( \mu \) et \( \sigma \).
La difficulté est le calcul de la vraissemblance marginale (en vert).
Il faut intégrer sur deux paramètres: \( \mu \) et \( \sigma \).
Remarque : la probabilité a posteriori est proportionnelle au produit de la vraissemblance et du prior.
Construction de la grille
n <- 200
d_grid <- tibble(
mu = seq(from = 150, to = 160, length.out = n),
sigma = seq(from = 6, to = 9, length.out = n)) %>%
expand(mu, sigma) # build all the combination -> matrix 2D
Fonction qui calcule la likelihood.
grid_function <- function(mu, sigma){
dnorm(d2$height, mean = mu, sd = sigma, log = T) %>%
sum()
}
On travaille en log pour réduire les erreurs d'approximation
d_grid <- d_grid %>%
mutate(log_likelihood = map2(mu, sigma, grid_function)) %>%
unnest() %>%
mutate(prior_mu = dnorm(mu, mean = 178, sd = 20, log = T),
prior_sigma = dunif(sigma, min = 0, max = 50, log = T)) %>%
mutate(product = log_likelihood + prior_mu + prior_sigma) %>%
mutate(probability = exp(product - max(product)))
head(d_grid)
# A tibble: 6 x 7
mu sigma log_likelihood prior_mu prior_sigma product probability
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 150 6 -1350. -4.89 -3.91 -1359. 1.91e-57
2 150 6.02 -1349. -4.89 -3.91 -1357. 5.73e-57
3 150 6.03 -1348. -4.89 -3.91 -1356. 1.69e-56
4 150 6.05 -1346. -4.89 -3.91 -1355. 4.95e-56
5 150 6.06 -1345. -4.89 -3.91 -1354. 1.43e-55
6 150 6.08 -1344. -4.89 -3.91 -1353. 4.07e-55
d_grid %>%
ggplot(aes(x = mu, y = sigma, z = probability)) +
geom_contour() +
labs(x = expression(mu),
y = expression(sigma)) +
coord_cartesian(xlim = range(d_grid$mu),
ylim = range(d_grid$sigma)) +
theme_bw(base_size = 30) +
theme(panel.grid = element_blank())
d_grid %>%
ggplot(aes(x = mu, y = sigma)) +
geom_raster(aes(fill = probability)) +
scale_fill_viridis_c(option = "A") +
labs(x = expression(mu),
y = expression(sigma)) +
coord_cartesian(xlim = range(d_grid$mu),
ylim = range(d_grid$sigma)) +
theme_bw(base_size = 30) +
theme(panel.grid = element_blank())
Cette méthode est beaucoup trop coûteuse et n'est pas utilisée en pratique
Pour échantillonner avec remise sur une grille, on utilise la fonction dplyr::sample_n()
set.seed(434)
d_grid %>%
sample_n(size = 1e4, replace = T, weight = probability) %>%
ggplot(aes(x = mu, y = sigma)) +
geom_point(size = 2., alpha = 1/15) +
scale_fill_viridis_c() +
labs(x = expression(mu[samples]),
y = expression(sigma[samples])) +
theme_bw(base_size = 30) +
theme(panel.grid = element_blank())
set.seed(434)
d_grid %>%
sample_n(size = 1e4, replace = T, weight = probability) %>%
select(mu, sigma) %>%
gather() %>%
ggplot(aes(x = value)) +
geom_density(fill = "steelblue", size = 0) +
scale_y_continuous(NULL, breaks = NULL) +
xlab(NULL) +
theme_bw(base_size = 30) +
theme(panel.grid = element_blank()) +
facet_wrap(~key, scales = "free")
library(tidybayes)
set.seed(434)
d_grid %>%
sample_n(size = 1e4, replace = T, weight = probability) %>%
select(mu, sigma) %>%
gather() %>%
group_by(key) %>%
mode_hdi(value)
# A tibble: 2 x 7
key value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 mu 155. 154. 155. 0.95 mode hdi
2 sigma 7.74 7.18 8.34 0.95 mode hdi
Le prior sur sigma n'est pas une distribution normale mais uniforme sur un intervalle.
Or on vient de voir que la distribution marginale postérieure semble normale
On peut se demander si c'est toujours le cas ???
Considérons le cas où le nombre de mesures est faible… ici 20!
set.seed(4341)
d <- read.csv("Data/Howell1.csv")
d2 <- d %>% filter(age >= 18)
d3 <- sample(d2$height, size = 20)
d3
[1] 161.925 147.955 156.845 161.290 163.830 161.290 153.035 159.385
[9] 139.700 154.305 146.050 154.940 149.860 152.400 153.035 160.655
[17] 161.290 148.590 144.780 154.305
1- définition de la grille
n <- 200
# note we've redefined the ranges of `mu` and `sigma`
d_grid <-
tibble(
mu = seq(from = 150, to = 170, length.out = n),
sigma = seq(from = 4, to = 20, length.out = n)) %>%
expand(mu, sigma)
2- définition de la fonction log_likelihood
grid_function <- function(mu, sigma){
dnorm(d3, mean = mu, sd = sigma, log = T) %>%
sum()
}
3- Calcul de la distribution postérieure
d_grid <- d_grid %>%
mutate(log_likelihood = map2_dbl(mu, sigma, grid_function)) %>%
mutate(prior_mu = dnorm(mu, mean = 178, sd = 20, log = T),
prior_sigma = dunif(sigma, min = 0, max = 50, log = T)) %>%
mutate(product = log_likelihood + prior_mu + prior_sigma) %>%
mutate(probability = exp(product - max(product)))
4- Échantillonnage sur la distribution postérieure
set.seed(4341)
d_grid_samples <- d_grid %>%
sample_n(size = 1e4, replace = T, weight = probability)
d_grid_samples %>%
ggplot(aes(x = mu, y = sigma)) +
geom_point(size = .9, alpha = 1/15) +
scale_fill_viridis_c() +
coord_cartesian(xlim = range(d_grid_samples$mu),
ylim = range(d_grid_samples$sigma)) +
labs(x = expression(mu[samples]),
y = expression(sigma[samples])) +
theme_bw(base_size = 30) +
theme(panel.grid = element_blank())
d_grid_samples %>%
select(mu, sigma) %>%
gather() %>%
ggplot(aes(x = value)) +
geom_density(colour = "steelblue", size = 1) +
scale_y_continuous(NULL, breaks = NULL) +
xlab(NULL) +
theme_bw(base_size = 30) +
theme(panel.grid = element_blank()) +
facet_wrap(~key, scales= "free")
La distribution sur sigma n'est plus vraiment Gaussienne…
Le modèle que l'on cherche à fitter sur les données
\[ \begin{align} h_{i} &\sim \mathrm{Normal}(\mu_,\sigma) \\ \mu &\sim \mathrm{Normal}(178, 20) \\ \sigma &\sim \mathrm{Uniform}(0, 50) \end{align} \]
Le code correspondant au modèle de gauche.
fit3.1 <- brm(
data = d2,
family = gaussian,
height ~ 1,
prior = c(
prior(normal(178, 20), class = Intercept),
prior(uniform(0, 50), class = sigma)
),
iter = 31000, warmup = 30000, chains = 4, cores = 4
)
La fonction qui fit le modèle sur les données est la fonction brm(…)
Le code qui permet d'afficher graphiquement le résultat
plot(fit3.1, theme = theme_bw(base_size = 30))
Le résultat de notre premier modèle en brms
plot(fit3.1, theme = theme_bw(base_size = 30))
Le résultat de notre premier modèle en brms.
Description complète
Family: gaussian
Links: mu = identity; sigma = identity
Formula: height ~ 1
Data: d2 (Number of observations: 352)
Samples: 4 chains, each with iter = 31000; warmup = 30000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 154.61 0.42 153.80 155.43 4165 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 7.76 0.29 7.21 8.34 817 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Description succinte baséée sur la moyenne et la variance (defaut)
posterior_summary(
fit3.1,
probs = c(0.025,0.975),
robust = FALSE,) %>%
round(digits = 2)
Estimate Est.Error Q2.5 Q97.5
b_Intercept 154.61 0.42 153.80 155.43
sigma 7.76 0.29 7.21 8.34
lp__ -1226.87 0.97 -1229.44 -1225.89
Description basée sur la médiane et la déviation absolue
posterior_summary(
fit3.1,
probs = c(0.025,0.975),
robust = TRUE,) %>%
round(digits = 2)
Estimate Est.Error Q2.5 Q97.5
b_Intercept 154.62 0.43 153.80 155.43
sigma 7.74 0.28 7.21 8.34
lp__ -1226.58 0.75 -1229.44 -1225.89
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 : Modèle linéaire à 1 prédicteur continu
Cours 4 : Modèle linéaire multivarié
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Prédicteurs catégoriels et Interactions
Cours 9: Modèles multi-niveaux
Cours 10: Data Hackaton